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ABSTRACT 

Cold accretion is a primary growth mechanism of simulated galaxies, yet observational evidence of 
“cold flows” at redshifts where they should be most efficient (z = 2-4) is scarce. In simulations, 
cold streams manifest as Lyman-limit absorption systems (LLSs) with low heavy-element abundances 
similar to those of the diffuse IGM. Here we report on an abundance survey of 17 H I-selected LLSs 
at z = 3.2-4.4 which exhibit no metal absorption in SDSS spectra. Using medium-resolution spectra 
obtained at Magellan, we derive ionization-corrected metallicities (or limits) with a Markov-Chain 
Monte Carlo sampling that accounts for the large uncertainty in N^i measurements typical of LLSs. 
The metal-poor LLS sample overlaps with the IGM in metallicity and is best described by a model 
where are drawn from the IGM chemical abundance distribution. These represent roughly 

half of all LLSs at these redshifts, suggesting that 28-40% of the general LLS population at 2 ^ 3.7 
could trace unprocessed gas. An ancillary sample of ten LLSs without any a priori metal-line selection 
is best fit with 48+of metallicities drawn from the IGM. We compare these results with regions of 
a moving-mesh simulation; the simulation finds only half as many baryons in IGM-metallicity LLSs, 
and most of these lie beyond the virial radius of the nearest galaxy halo. A statistically significant 
fraction of all LLSs have low metallicity and therefore represent candidates for accreting gas; large- 
volume simulations can establish what fraction of these candidates actually lie near galaxies and the 
observational prospects for detecting the presumed hosts in emission. 

Subject headings: galaxies: evolution, intergalactic medium, high-redshift, quasars: absorption lines 


1. INTRODUCTION 

Over the past decade, a new theoretical paradigm de¬ 
scribing galaxy evolution and gas accretion has emerged 
from the synergy between semi-analytic galaxy forma¬ 
tion modeling and high-redshift observations. In pre¬ 
vailing models of galaxy formation, spiral galaxies grow 
largely through a hierarchical merger process, with rel¬ 
atively quiescent star formation driven by gas accretion 
onto the dark matter halo and major mergers initiat¬ 
ing periods of rapid starburst, ultimately resulting in 
elliptical galaxies with quenched star formation (Kauff- 
man n et al.||1993l |Mo et al.||1998| |Toomre & Toomre 
1972[ Hopkins et al.|2007). However, recent morphologi¬ 
cal evidence indicates that disk galaxies at high redshifts 
grow largely through smoo th gas accretio n di rect ly onto 
the stellar disk ( Bonrnand fe Elmegreen|2QQ9 Bournaud 
et al. 2009), and mergers may play a less prominent role 
in their growth (van Dokkum et al. 2013). Additionally, 
the major-merger rate is too low to explain the num¬ 
ber of gal axies at z ~ 2-3 with a high star-form ation 
rate (SFR,|Jogee et al.|2008||Elmegreen et al.|2007|), and 
such mergers do not predict m orphologies seen m galax¬ 
ies with high SFRs (Dekel et al.||2009a). Furthermore, 
a class of massive compact spheroidal galax ies with low 
SFRs is already well established at z ~ 2 (Kriek et al 
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van Dokkum et al. 2008), inconsistent with a sce¬ 


nario in which they are solely a product of major mergers 
between disk-like galaxies. 

Recent simulation work has explored a complemen¬ 
tary galaxy growth mechanism, in which massive galax¬ 
ies at high redshift are stream-fed large quantities of gas, 
and properties of the accreting g as influence the resul - 
tant SFR and morphology (e.g., |Dekel et al.||2009a|b|). 
Central to this framework is the existence of “cold-flow” 
accretion—filamentary gas traveling directly from the in¬ 
tergalactic medium (IGM) onto the star-forming disks of 
galaxies, without shock heating at the virial radius. Al- 


accreting via cold flows (e.g., 

Nelson et al. 

2013 

), it re- 

mains a common feature, fundamental to t 

he growth of 


early star-forming galaxies. 

Observationally, cold flows are expected to manifest as 
Lyman-limit systems (LLSs), absorption systems along 
quasar sightlines with 7912 > 2 (logTVni > 17.5), in the 
extended intra-halo medium of galaxies. This environ¬ 
ment, often referred to as the circumgalactic medium 
(CGM), forms a regulatory interface between galaxies 
and the IGM and p otenti ally holds a large reservoir of 


baryons (Werk et al. 


2013). 


Although the CGM also contains outflowing and recy¬ 
cling gas from the galaxy t hat manifests as absorption 
(e.g., Bordoloi et al. 2014), the gas metallicity of ab¬ 
sorbers serves as a straightforward diagnostic to distin¬ 
guish between accreting baryons and other gas. Out¬ 
flowing and recycli ng gas have been enriched to high 
metallicities (|Steidel et al. ||2010), often approaching so¬ 
lar, whereas gas that isoemg newly introduced to a 
galaxy is more likely to have low metallicities consistent 
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with that of the diffuse IGM (Fumagalli et al. 2011b). 
As such, observations establishing the prevalence of LLSs 
with pristine elemental abundances would provide signif¬ 
icant evidence supporting cold-flow models of baryonic 
accretion. 

Cold-flow accretion should be most efficient at 2.5 < 


z < 4. 5, during the peak of cosmic star formation (Keres 
et al. 2009); indeed simulations find significantly larger 
cover i ng fractions at z ^ 4 (Faucher-Giguere & Keres 
2011 Kimm et al. 2011). Several groups have reported 
detections of individual low metallicity LLSs indicative 
of cold-flow accretion, but the larg er population of high-z 


LLSs remains mostly unexplored. Crighton et al. (2013b 
report the discovery of a z = 2.44 LLS with metallicitw] 
[M/H] = —2.00±0.17 nea r a low-luminosity galax y mixed 
with metal-rich material; |Levshakov et al. (20 03|) discuss 
a z = 2.92 LL S with [C/HJ = —2.93±0.13; and|Fumagalli 
(|2011a) present two LLSs at z = 3.41 (3.10) with 
limits ( 


et al 


provides a large, H I-selected high-redshift sample that 
can be used to perform such a study on the population 
of high-redshift LLSs. 

For this paper, we constructed a survey of high- 
redshift, H I-selected LLSs exhibiting low heavy-element 
abundances in SDSS spectra. In Section 2, we discuss the 
selection of sight lines, new observations, and data pro¬ 
cessing. Our measurements, ionization modeling, and 
metallicity determination of the LLSs are described in 
Section 3. Section 4 presents the measured metallicities 
of the observed LLSs. In Section 5, we discuss prop¬ 
erties of the LLS population, quantify the fraction of 
low-metallicity LLSs that are candidates for the obser¬ 
vational signature of cold-flow accretion, and compare 
with simulations and other observational studies. Section 
6 provides a summary. Throughout, we adopt ACDM 
cosmology with cosmo logical parameters from WMAP9 
(Hinshaw et al. 2013): = 0.72, = 0.28, and 


upper 


of [M /H] < -4.2 (-3.8). Additionally, two Hq = 70 km s 1 Mpc 


of thes e absorbers (|Crighton et al.||2013 Fumagalli et al. 
2011a) have clear deuterium detections with column den- 


sities consistent with primordial abundances, indicating 
the gas comprising the abs orbers has had little m ixing 
with gas processed by stars. ]]Fumagalli et al. (2013) con¬ 
struct a composite absorption spectrum from 20 LLSs at 
z ~2.6-3 selected from a blind QSO survey and find, for 
the composite, [M/H] < —1.5, similar to the observed 
metallicities of damped Ly<a systems (DLAs, absorbers 
with logTVni > 20.3). 

Studies of the LLS p opulati on and CGM are more ex¬ 
tensive at low redshift. Lehner et al. (2013) study 28 H I- 
selected LLSs at z < 1 and find a" bimodality in metal- 
licity, with peaks at [M/H] ~ —1.60 and —0.3. Addi¬ 
tionally, studies connecting low-redshift metal absorbers 
with host galaxies find that the distribution of absorbers 
is azimuthally and morphologically dependent in a fash¬ 
ion consistent with a general picture of ga s accretion and 
galactic winds (e.g., Kacprzak et al. 2012), although such 
studies have not yet been coupled to metallicity. The 
low-metallicity branch of the bimodal distribution is con¬ 
sistent with the notion of cold gas reservoirs. However, 
star formation and accretion rates are much lower during 
this epoch than at highe r redsh ifts. 


In simulations, Neistein et al. (2006) found the average 
accretion rate onto galacti c halos m ACDM c osmology 
goes as M oc (1 ± z) 2,25 . Dekel et al. (2009b) showed 
that baryonic-input rates from cold gas streams in cos¬ 
mological hydrodynamical simulat ions follow the same 
expression at high redshift. If the Lehner et al. (2013) 
metallicity bimodality reflects a distinction between ac¬ 
creting gas and enriched outflowing or recycling gas, then 
it should be more pronounced during the peak of cosmic 
star formation. 

There is claimed evidence of cold-flow accretion at low 
redshift, and theoretical predictions indicate an increas¬ 
ing frequency with redshift. However, there are only a 
handful of high-redshift observations indicative of cold 
flows. We seek to determine whether this is due to an 
observational shortage or a departure from expectations 
from simulations. The Sloan Digital Sky Survey (SDSS) 

5 Metallicity is denoted as [X/H] = logfVxAVn) — 

log(iVx,© Ah,©)) where iVx is the column density of an arbitrary 
atomic species. Often we report [M/H] to indicate “all metals.” 


2. DATA 


We pre-screened a large, H I-selected sample of LLSs 
for candidates likely to have low metallicities and ob¬ 
tained higher resolution spectra of 15 candidates. This 
approach maximized information about the metal-poor 
end of the LLS distribution, which was not well under¬ 
stood at the outset. However, it also split the statisti¬ 
cally characterized H I sample, complicating the broader 
interpretation of results. The pre-screening proved less 
efficient than expected in identifying ideal candidates, 
roughly halving the parent sample. We focused our 
first observations on this metal-poor sub-population, and 
work exploring the full LLS population is underway. 

The parent sample is 194 LL Ss with zt.t.s > 3 .3 and 
Nui > 17.5 cm 2 , compiled by |Prochaska et al. | (|2010| 
hereafter POW10) using quasar spectra from the SDSS 
Data-Release 7 (DR7). They identified systems by con¬ 
structing models of the Lyman-series absorption and the 
Lyman break and applying them to absorbed quasar 
continua. Although they identify more than 194 LLSs, 
we only consider their “statistical sample”, comprised of 
spectra with zqso > 3.6 and zlls > 3.3. 

We cross-referenced this list with a C IV AA1548,1550 
catalog co nstruc ted from the SDSS DR7 quasar spectra 
([Cooksey et al.|2013|) and found that 152 of the 188 SDSS 
spectra in our parent sample did not have associated C IV 
detections within ±500 km s -1 of the LLS redshift. We 
visually inspected the remaining 152 spectra for typi¬ 
cal metal absorption lines at zlls- Several spectra had 
weak C IV doublets below thresholds of the C IV sur¬ 
vey or C IV obstructed by interlopers, and many spectra 
did not have definitive C IV but had absorption from 
C II A1334 or other low-ionization species. Ultimately, 
the LLS sample was roughly halved by a metal-line in¬ 
spection, with definite or probable metal absorption lines 
associated with 100 of the LLSs, and no corresponding 
metals seen for 96 of them. 

The metal-poor LLS sample presented here was sub¬ 
jected to an additional declination cut for observation 
at the Magellan telescopes, since they are situated at 
a latitude of -29° and the SDSS footprint is primarily 
in the northern sky. Excluding quasars with S > ±21° 
(corresponding to a transit airmass of « 1.5) leaves 28 
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sight lines, 15 of which we observed in this initial survey. 

2.1. MagE Spectra 

We obtained higher resolution spectra along 15 quasar 
sightlines (Tabled selected as described ab ove usin g the 


Mage llan Echellete Spectrograph (MagE, Marshall et al. 
2008) on the 6.5-m Magellan Clay telescope. MagE cov- 


ers optical wavelengths from 3100 A to 1 /im. At zlls 
3.3, the 912 A Lyman break is redshifted to 3900 A. With 
an 0.85" slit, MagE has a resolution 7 Z = 4950 (or full 
width at half maximum FWHM = 60.7km s -1 ). Obser¬ 
vations were done on UT 17/19 March 2013 and UT 05 
May 2013 with typical seeing of 0.6 // -0.8 // . 


Data were reduced using the MASE pipeline (Bochan- 
ski et al. 2009), using GJ 620.1 B/HIP 80300 as a stan¬ 


dard for calibration. MASE is an IDL software package 
designed for reducing MagE data and performs the full 
extraction and calibration process. We manually con¬ 
struct a cubic-spline fit to the continuum of each spec¬ 
trum. 

Figure [l] shows portions of the SDSS and MagE spec¬ 
tra of J083832 around several of the LLS metal lines for 
comparison. In this example there is no statistically sig¬ 
nificant metal absorption from the LLS in the SDSS spec¬ 
trum (FWHM ~ 150 km s -1 ) but the MagE spectrum 
has clear absorption lines. 

2.2. Higher Resolution and Infrared Spectra 

Several of the LLSs have no metal absorption in their 
corresponding MagE spectra. We observed one such ob¬ 
ject, J124957, with the Ma gellan Inamori Kyocera Echel- 
lette spectrograph (MIKE, |Bernstein et al.||2003|) at the 
same telescope on UT 20 March 2013 using a 1" slit. 
MIKE is a double echelle spectrograph. The blue arm 
covers 3200 A to 5100 A, and the red arm covers 4900 A to 
1 /rm. With a 1" slit, the blue/red arms have resolutions 
of 28,000/22,000 (FWHM = 10.7kms -1 / 13 - 6 kms -1 ). 
All metal lines used in the current survey are in the red 
portion of the spectrum. MIKE data were reduced using 
MIKE ReduxJ^] a series of IDL tools that encompass all 
calibrations and extractions. 

We also make use of a medium-resolution infrared spec¬ 
trum of the same object taken with t he Fol ded-port In- 
fraRed Echellette (FIRE,|Simcoe et al.| 2008), on the 6.5- 
m Magellan Baade telescope, observed as part of a differ¬ 
ent survey. FIRE has a bandpass covering 0.8/im-2.5/im, 
at a resolution of 1Z = 6000 (FWHM = 50km_s_ -1 ). 
For data acquisi tion and reduction details see Matejek 

& Sim^(|20T2). 

The primary motivation for the MIKE observation was 
the possibility of identifying deuterium absorption asso¬ 
ciated with the H I, wh ich can indicate that gas is unpro¬ 
cessed (e.g., Crighton et al. 2013). Unfortunately, the H I 
absorption from the LLS along this sightline proved too 
broad to distinguish deuterium absorption. However, the 
higher resolution MIKE spectrum enables more sensitive 
column-density measurements, and the FIRE spectrum 
allows us to measure several ions not covered by the op¬ 
tical instruments. 


2.3. Metal-Blind Sample 


6 http://web.mit.edu/~burles/www/ 


SilV 1393 Sill 1526 CIV 1548 CIV 1550 


nJU.—n-uj- 

11—1 — i-H n 1—1 —■ —i 

n—T1 r—i r—i n 


SDSS 

L-l 1 1 j| 

l— n i-h-lJ 


''' i'' 1 1 1111111 

MagE 

''i| rrn~rrrprrr 

"Jl i ii 1 ii i h i i 

- 


-500-250 0 250-500-250 0 250-500-250 0 250-500-250 0 250 

Relative Velocity (km s' 1 ) 


Fig. 1.— Comparison of normalized spectra (black) of J083832 
from SDSS and MagE. The red line is the la error; the green line is 
the continuum (unity). The spectrum cutouts are centered around 
^lls (v = Okms -1 ). The Si IV A1393 and Si II A1526 lines and 
C IV doublet are all evident in the MagE spectrum, but are not 
seen in the SDSS spectrum. Si IV A1402 is unavailable in both 
spectra due to a strong interloping absorber. 


We are interested in how our sample of metal-poor 
LLSs compares to the global LLS population. As a com¬ 
parison s ample, we studied spect ra from the blind LLS 
survey of Fumagalli et al. (2013), also conducted with 
MagE. Since their sample is at slightly lower redshifts 
than ours, we selected the ten highest-redshift absorbers 
from their survey (excluding several DLAs and absorbers 
close to the QSO redshift) to achieve the best redshift 
overlap with our metal-poor sample. This method of 
choosing objects also avoids introducing any selection 
bias with regards to metallicity. The ten LLSs we ex¬ 
amined have a median redshift of £lls=3.04. For the 
remainder of the paper, we refer to this dataset as the 
“metal-blind sample,” and to our observations as the 
“metal-poor sample.” We applied identical analysis tech¬ 
niques to reduced spectra in both samples. 


3. ANALYSIS 

We inspected the Lyman series absorption in each of 
the MagE spectra to confirm LLS redshifts found in 
POW1Q using SDSS spectra. We found a difference in 
redshift of <0.01 for all but one of the systems. The 
outlier J085944 has a weaker absor ber (logTVni < 17.5) 
at the redshift determined by |PQW10[ The redshift of 
the LLS (^lls = 3.263) is smaller by « 0.2 and is the 
lowest redshift LLS in our study. Since we did not check 
for metal absorption in the SDSS spectrum at the cor¬ 
rect redshift, this system could have biased our sample. 
However, the absorber serendipitously met the target se¬ 
lection criteria discussed in Section |2j so we included it 
in our analysis. Additionally, two oft he quasar spectra 
had a second, lower-redshift LLS close enough in redshift 
to target system that enough of the Lyman series transi¬ 
tions were available to measure the H I column density. 
These two LLSs also meet the selection criteria (no met¬ 
als seen in the corresponding SDSS spectrum) and were 
included in the analysis. 


3.1. HI Column Density 

TVhi is notoriously difficult to measure in the LLS 
regime; the Ly<a curve-of-growth is flat and higher order 
transitions, including the Lyman break, are saturated by 
construction. We found that for the purpose of calcu¬ 
lating metallicities, ionization modeling techniques can 
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H0449 (17.8 <N h ,< 18.2) 



J093153-000051 (18.4 < N HI < 19.3) 



J085944+212511 (18.8 < N H , < 19.6) 
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Fig. 2.— Velocity plots of example H I profile fits. The black curves are the normalized spectra, and the red curves are the la error on 
the flux. Green lines are representative fits, and the cyan shadings fill the regions between the Voigt profiles corresponding to the smallest 
and largest plausible H I column densities. Left: Several Lya absorption profiles, showing how Lya fitting contributes in different regimes 
of Nm. The bottom profile has damping wings that tightly constrain H I. The figures above this show how Lyo; can place an upper limit 
on Nm by requiring model profiles to not over-absorb. Center: Several Lyman series transitions drawn from the same system. For this 
absorber, Higher order Lyman series transitions (e.g., Ly6 A930, Ly8 A923, Lyl2 A917) are used to measure the Doppler b parameter, and 
lines nearing the Lyman limit are fit to measure Nm (top panel). The saturation at the Lyman limit also places a lower limit on JVhi. 
Right: Example of a comparatively weak absorption system. The absorption is not fully saturated at the Lyman limit, allowing an Nm 
estimation. The absorption profile does not vary within the range of Nm allowed by the flux seen at the Lyman limit (top panel). 


marginalize over a wide range of H I column densities, 
at least within the LLS range, as we discuss in Section 

m 

Rather than attempting to find H I explicitly using 
Voigt profiles, we determined a range of viable column 
densities for each system, listed in Table [2] We used mod¬ 
ified versions of IDL software from the XIDlj^] library in 
conjunction with the Voigt profile fitting packages VPFIT 
and RDGEI0 

For each system, the plausible range of Nm is deter¬ 
mined by fitting various different aspects of the H I ab¬ 
sorption signature (Figure [2]). We estimated the Doppler 
b parameter, a measure of the width of the Voigt profile, 
using higher order Lyman series transitions. For weaker 
systems where the Lyman limit is not fully saturated, 
we were able to constrain Nm using the flux at the Ly¬ 
man limit and/or fitting higher order Lyman series lines. 
For the strongest absorption systems, we fit the weak 
damping wings on the Ly<a profile. For systems of mid¬ 
dling strength, the saturation of the Lyman limit coupled 
with the non-existence of Ly<a damping wings restrict the 

7 http://www.ucolick.ors/~xavier/IDL/ 

8 http://www.ast.cam.ac.uk/~rfc/vpfit.html 


range of possible Nm. 

The typical Nm uncertainty (defined as max(VHi)- 
min(TVHi) for a given system) for both samples is 0.7 
dex. This median total deviation is akin to an error bar 
of ±0.35 dex. The best constrained system had a total 
deviation of 0.3 dex, while the least constrained had a 
deviation of 1.7 dex, although the maximum is an outlier 
with the next largest being 1.3 dex. These errors are 
incorporated into our metallicity uncertainty as bounds 
on a flat prior distribution of allowed Nm , defining the 
range where we explore possible values for our solutions. 


3.2. Metal Column Densities 

We measure column densities for ionic spe cies using 
the apparent o ptical depth method (AODM, | Savage fc 


Sembach 1991). For each absorber, we integrated over 


a fixed velocity width in order to maintain consistency 
between different species/lines (rounded to the nearest 
pixel). For ions with multiple available lines, we per¬ 
formed an error-weighted average of the detections. 

For each line where there was a non-detection, we cal¬ 
culated a 3 -<7 upper limit to the column density corre¬ 
sponding to the error in the spectrum over ±1 resolu- 
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tion element using 10,0000 Monte-Carlo realizations. For 
each iteration, we added to the flux a value drawn from 
a Gaussian distribution with width equal to the error 
at that pixel, then measured the column density of the 
resulting mock profile, discarding realizations where the 
column density was negative. If an equal amount of flux 
were scattered above and below the continuum level, this 
would result in a positive column density because the 
relationship bet ween flux and apparent optical depth is 
exponential (see Fox et al.|2005 ). This produced a distri¬ 
bution of the largest possible column densities consistent 
with the observed lack of absorption. We adopted the 
column density larger than 99.7% of other trials as the 
3 -a upper limit. To mitigate the effects of poor contin¬ 
uum fits resembling low-column density absorption, we 
set the flux to unity and repeated the process if less than 
50% of the trials produce a positive column density. 

Similar to the H I analysis, we dismissed all metal ab¬ 
sorption lines contaminated by interloping absorbers as 
well as lines obscured by large amounts of noise. The 
Si II A1260 transition, which is a powerful diagnostic 
for low-metallicity absorbers due to its large oscillator 
strength, was unavailable for ten of the 17 systems due 
to confusion with the Ly<a forest. 

All AOD column density measurements assumed the 
absorption profiles reside on the linear portion of the 
curve-of-growth and are unsaturated. While this as¬ 
sumption is expected to hold for all of the lines in our 
metal-poor sample, there is no absorption strength cut 
for the metal-blind sample so we need to test for satura¬ 
tion. In low- and medium-resolution spectra, absorption 
profiles can be saturated without clearly exhibiting zero 
flux, since the instrument blurs the absorption profile. 
In Table [l] we list the measured column densities, and 
below we discuss the identification of saturated lines. 


3.2.1. Testing for Saturation 

The AODM provides a test for saturation through 
comparison of the the AOD profi les fo r different tran¬ 
sitions of the same species (Savage & Sembach 1991). 
However, this is insufficient when a species only has a sin¬ 
gle, potentially saturated line. To test for saturation in 
such species, we performed multi-component Voigt pro¬ 
file fits to see if we recover column densities similar to 
those measured with AODM. Since the velocity structure 
of absorbers is typically not well resolved in our spectra, 
we used Monte Carlo methods to explore the parameter 
space for each line. For each line, we constructed 200 
models to use as input to VPFIT, each having 3-5 com¬ 
ponents (uniformly selected) placed by splitting the ab¬ 
sorber into bins and uniformly selecting redshifts within 
these bins, such that components extend over the en¬ 
tire absorption profile. Foll owing measurements of C IV 
Doppler b parameters (Rauch et al. 1996), components 
have Doppler parameters drawn from a Gaussian dis¬ 
tribution with b = 12kms _1 and oq = 5kms _1 , con¬ 
strained to be above 6kms _1 , with b held fixed during 
fitting. We allowed VPFIT to remove components and 
slightly modify redshifts. 

The intent of this exercise was not to determine spe¬ 
cific models for the unresolved velocity structure of the 
absorbers but rather to gauge whether there is satura¬ 
tion by seeing if the models with column densities larger 


TABLE 1 

Details on specific absorbers 


Ion Arest logAAoDM log ^adopted"* 


Metal-Poor Sample 


J080853-070940 s = 3.545 


C IV 

1548 

13.38 ±0.13 

13.39 ±0.11 

C IV 

1550 

13.42 ±0.21 


O I 

1302 

< 13.65 

< 13.65 

Si II 

1304 

< 13.31 

< 13.30 

Si II 

1526 

< 13.45 


Si IV 

1393 

12.93 ±0.12 

12.93 ±0.12 

Si IV 

1402 

< 13.36 



J083832+200142 s 

= 3.47595 

Al II 

1670 

12.73 ± 0.04 

12.73 ±0.04 

C II 

1334 

14.29 ± 0.02 

> 14.29 

C IV 

1548 

13.65 ±0.05 

13.70 ±0.04 

C IV 

1550 

13.85 ±0.07 


Si II 

1526 

13.49 ±0.10 

13.49 ±0.10 

Si IV 

1393 

13.52 ±0.03 

13.52 ±0.03 


Note. — This table will be published in 
its entirety in the electronic edition; a portion 
is shown here as an example. 


than measured by the AODM provided reasonable fits 
to the spectra. In evaluating the output, we first re¬ 
moved all fits where the structure was reduced to a sin¬ 
gle component (which tends to produce poor fits with 
unrealistically large column densities, since the Doppler 
parameters were fixed) and all fits where the x 2 -fit statis¬ 
tic output by VPFIT was more than 2 a above the mean 
(these are typically trials where the profile is essentially a 
single component fit with several negligible components). 
The remaining models provide a distribution of potential 
column densities for the absorbing ion. 

We considered an ion to be saturated if it met two re¬ 
quirements: (i) less than 5% of the trials had column 
densities less than that obtained via AODM and (ii) the 
median column density of the modeled distribution ex¬ 
ceeded Vaodm + max(3<TAODM 5 0.2 dex). The second cri¬ 
terion prevented lines with Monte Carlo trials resulting 
in very precise, narrow output distributions only slightly 
larger than the AODM measurement from being falsely 
classified as saturated. When we determined that an ion 
was saturated, we adopted the AODM measurement as 
a lower limit to its column density. 

We found this approach agreed with both our expecta¬ 
tions for which lines are saturated based on appearance 
and AODM testing for saturation. For very weak ab¬ 
sorption profiles this technique did not reliably produce 
meaningful results, but comparing AOD profiles when 
there are multiple ions shows that such lines are clearly 
unsaturated, as expected. In our metal-blind sample, we 
found only one saturated line, a C II A1334 line that was 
misattributed to QSO H I self-absorption in the SDSS 
spectrum, both due to its location on the QSO Ly<a emis¬ 
sion peak and the lack of corroborating lines. There are 
numerous saturated lines in the metal-blind sample, as 
indicated in Table [lj 

3.3. Ionization Modeling and Metallicity 
Determination 
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Fig. 3.— Absorption profiles and example ionization model for 
J123525. Left: The normalized spectrum (black) of J123525 at 
the position of various metal lines relative to zlls ( v — Okms -1 ), 
with lcr uncertainty on the flux in red. The green line indicates the 
continuum (unity). The blue lines indicate the range over which 
the optical depth was integrated to determine the column density. 
Where no absorption was detected, the blue lines are instead d= one 
resolution element from the central redshift, indicating the range 
over which the 3a upper limit was measured. Right: Example 
ionization model. For each ion, the solid and dashed curves show 
the column density as a function of the ionization parameter for 
different metallicities. The blue shaded regions are lcr intervals 
around the column density detections or, for non-detections, the 
region below the upper limit. The yellow shaded region is the range 
of log U found to have a consistent solution. The black line is the 
ionization parameter corresponding to the best-fit solution. 


The primary interest of our study is the metallicity, 
[M/H], of the absorbers. However, this requires knowl¬ 
edge of elemental abundance ratios, rather than the in¬ 
formation on specific ions that we measure. Rather than 
assuming ionization conditions to convert between ionic 
and atomic column densities, we used the softwa re pack¬ 
age Cloudy (version 13.02, last described by, Ferland 


2013b to solve for the ionization and metallicity 


et al. 

simultaneously. With Cloudy, we modeled the ioniza¬ 
tion conditions of the absorbers over a range of metal- 
licities, obtained ionic column densities for each model, 
and determined which models best matched the observed 
column densities using Monte Carlo simulations. Ioniza¬ 
tion conditions are typically described by the ionization 
parameter [/, a proxy for hydrogen density n h defined 
as: 


U= — 
n H c 


(i) 


where c is the speed of light and nn = i + ii + ^h 2 • 
The flux of H I-ionizing photons, <L, is given by 


T = 4tt 



^dv 

hv 


( 2 ) 


Fig. 4.— Same as Figure [3] but for J124957. Left: The normal¬ 
ized spectrum of J124957 at the position where it’s LLS’s various 
metal lines would be. Since no absorption is detected for any of 
these lines we only obtain upper limits. Right: Without any metal 
column density measurements to constrain the ionization param¬ 
eter, we assume log U > —3 and measure the metallicity upper 
limit at log U = — 3. For each ion, the limiting metallicity is at the 
upper-left corner of the overlapping shaded regions; the lowest of 
these gives the upper limit metallicity for the LLS. For this system, 
the strongest constraint at log U = —3 comes from Si II 1260. 


where J u is the specific intensity of the incident radiation 
and z^ll is the frequency corresponding to 1 Ryd. 

Specifically, we used Cloudy to calculate the column 
densities of different species as a function of log V over a 
grid of H I column densities, redshifts, and metallicities. 
The models were generated assuming a plane-parallel ge¬ 
ometry of a uniform and isothermal layer of photoionized 
gas, with the shape of the ionizing radiation spectrum de¬ 
rived from a combination of the cosmic microwave back¬ 
ground (CMB) and th e cosmic ultraviolet bac kground 
(UVB) spectrum from |Haardt & Madau| (|2012| CUBA 
software). The CMB is much weaker than the UVB at 
all relevant wavelengths and omitting it caused no ap¬ 
preciable change. We adopted a solar relative abundance 
pattern for the Cloudy models. 

The UVB spectrum includes emission from galaxies 
and QSOs, as well as a saw tooth absorption pattern due 
to the He II Lyman series ( Madau fc Haardt| [2009 ). Al¬ 
though observations of the hydrogen ionizing fl ux dis 


agree with the normalization of this spectrum (Becker 


& Bolton 

2013 

Faucher-Giguere et al. 2008), Cloudy 

models ad 
value of V 
the shape 

just the normalization according to the input 
r , so this is inconsequential. Efforts to adjust 
of the spectrum have found that best-fit mod- 


els typically require fairly small modifications that ulti- 
mate ly translate to a difference in metallicity of < 0.2 
dex (Crighton et al. 2015), although there are outliers 
that are best fit with larger modifications to the shape. 
Since our modeling is based on a small number of ions, we 































































































































































Low-Metallicity LLSs 


7 


keep the number of fit parameters to a minimum to avoid 
overfitting and do not let the shape of the spectrum vary 
in our models. Additionally, this ionizing background 
(without adjustment) was u sed in most of the works we 
compare to in Section |5.3[ so any uncertainties in the 
shape of the spectrum should minimally affect compari¬ 
son with other observations. 

For each LLS, we interpolated on this grid to the ab¬ 
sorber redshift, producing for each ionic metal species 
the column density as a function of [M/H], log [/, and 
TVhi. In Figures [3] and [4] we sketch the ionization model¬ 
ing technique for LLSs with and without metal-line ab¬ 
sorption, respectively. These examples assume a value 
of TVhi intermediate to the allowed range, so one dimen¬ 
sion is missing from this schematic representation of the 
modeling process. 

3.3.1. Dependence of Results on TVhi 

Although many systems require us to marginalize over 
a fairly broad range in TVhi, ionization modeling within 
the bounds we consider is surprisingly insensitive to the 
particular H I column density. Just as Figures [ 3 ] and [4] 
are projections onto an assumed value of TVhi, in Figure 
[ 5 ] we project along two different values of log U (at fixed 
[M/H] = —2.5), to clarify how several properties vary 
with TVhi- 

In the top panel, we plot the model column density as 
a function of TVhi for several different ions, as well as the 
total hydrogen column density TVh, scaled to fit on the 
same plot. For a given value of log [/, TVh does not scale 
very rapidly with TVhi; over the two orders-of-magnitude 
of TVhi shown, TVh only changes by about 0.3 dex. The 
column density curves for C IV and Si IV are comparably 
flat. The variation in the low-ionization metals with TVhi 
is more appreciable, but is still one order-of-magnitude 
smaller than the variation in TVhi- We do not show the 
scaling with [M/H], although it is as expected-the metal 
column densities increase by an order-of-magnitude when 
[M/H] is increased by one. Since TVh and most of the 
metal column density curves are relatively flat functions 
of TVhi, the model parameters V and [M/H] correspond¬ 
ing to the model that matches measured ionic column 
densities from an LLS do not vary strongly with TVhi- 

Comparing the model column density curves for the 
two different ionization parameters plotted, we see that 
the low ions are not strongly dependent on log U. C IV 
and Si IV column densities depend much more strongly 
on how ionized the gas is, as can also be seen in the 
right-hand plots of Figures [ 3 ] and [4] 

We estimate from the Cloudy output how a metallicity 
measurement (at fixed log U) based on a single ion varies 
with TVhi- For ion x corresponding to atom X, with mea¬ 
sured column density N x and model ionization fraction 
f x = TVe/TVx, we can write the metallicity as 

[X/H] = log TV*/TV h - log(TV x /TV H )© 

= log - log(Nx/N H ) Q 

Ac tv h i 


where /hi is the H I ionization fraction in the correspond¬ 
ing model. Noting that the measured column density is 



17.5 18.0 18.5 19.0 19.5 

log N m (cm' 2 ) 


Fig. 5.— Top: Model column densities for several ions, assum¬ 
ing logV = — 2.2 (solid curves) and log U = — 2.6 (dashed curves) 
with fixed metallicity [M/H] = —2.5. The black curve is the total 
hydrogen column density TV#, scaled by 10 -8 . The intermediate 
ions (C IV and Si IV), along with TVh are quite flat with TVhi- 
The low ions are more correlated with TVh I although the range of 
column densities they span is still about one-tenth of that consid¬ 
ered in TVhi- Bottom: Slope of the metallicity as measured by a 
single ion and assuming a fixed value of log U. This can be used to 
estimate the uncertainty in [M/H] introduced by the uncertainty 
in TVhi, although this approach overestimates [M/H] uncertainty 
since it does not consider overlapping constrains from multiple ions. 

a constant, we can differentiate with respect to log TVhi: 
0[M/H] = 8 log / m d log/s t 
d log TV h 1 9 log TV H 1 d log TV H 1 

The derivative depends only on ionization model output, 
and is shown in the bottom panel of Figure [3] for several 
ions. 

We first note that all of the slopes are negative, indicat¬ 
ing that metallicity decreases with TVh i as expected. The 
Si IV and C IV slopes approach zero above TVhi ^ 18, in¬ 
dicating that modeled metallicities should be very robust 
in this range of column densities—given TVqiv or TVsnv 
and log£/, the metallicity is independent of TVhi- Si II 
and C II have larger derivatives with TVhi, as expected 
from the column density curves, but weaker dependence 
on log U. 

This suggests that uncertainty of TVhi is manageable. 
Treating the metallicity slope as a proxy for the model 
uncertainty in the resulting metallicity and integrating 
the C II derivative from TVhi=17.5 to 18.5, where the 
model uncertainty is worst, results in only a ^ 0.5 dex 
uncertainty in metallicity. In practice, the uncertainties 
are significantly less since models are fit using multiple 
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ions. 


3.3.2. Matching Ionization Models to LLSs 

To compare how well different ionization models fit the 
data, we define the likelihood function 

£([M/H], TVhi, logt/) = Hi£i(Ni) 

where £i(Ni) is the likelihood for each ion measured given 
N i: the model column density obtained by interpolating 
the grid to the corresponding values of [M/H], TVhi and 
log U. For detections, li is taken to be a Gaussian with a 
mean and standard error given by the AODM measure¬ 
ments. For upper- and lower-limits, we let li be unity 
if the model column density is within the range allowed 
by the limits, and li decays as a Gaussian with a = 0.05 
beyond the allowed range. In practice, we used log£ 
to avoid computational instabilities resulting from small 
likelihoods. _ 

Motivated by the implementation of [Crighton et al 
(|2015|), we used th e Python module emcee (|Foreman- 
Mackey et al. 2013) to perform a Markov Chain Monte 
Carlo (lVlCMC) sampling of the parameter space. This 
approach allows for calculation of the posterior proba¬ 
bilities of the metallicity and ionization parameter while 
marginalizing over the possible values of TVhi- We used 
flat priors with log U, [M/H] E [—4, —1], and TVhi within 
the viable range determined for each LLS. We ran 1000 
iterations of 100 walkers sampling the parameter space, 
discarding the first 100 iterations as ‘burn-in’ to allow the 
walkers to explore the full posterior distribution and to 
remove the signature of the walkers’ initial conditions. 
We constructed the posterior distribution from the re¬ 
maining 900 iterations. An example posterior distribu¬ 
tion (for J123525) is shown in the top left of Figure [6j 

Although many systems required marginalization over 
a wide range in TVhi, often larger than one order-of- 
magnitude, we found that the ionic column densities 
and posterior probabilities are not strongly dependent on 
TVhi, c onsiste nt with expectations from the discussion in 
Section l3.3.11 


3.3.3. Metallicity Upper Limits 

For LLSs where the absorption line data were insuffi¬ 
cient to constrain the solution, we derived upper limits 
to the metallicity. The way an upper limit was found 
depends on whether or not there were any metal-line de¬ 
tections. We note that as log U increases, nn decreases 
and the gas becomes more highly ionized, resulting in col¬ 
umn densities for the ions we measure to correspond to 
lower metallicities at larger ionization parameters (Fig¬ 
ures [3] and |4j and Figure [6] left-middle panel). Hence, 
lower values of log U correspond to more conservative 
upper bounds on metallicity. 

For absorbers with only non-detections, we found the 
limiting metallicity at log U = —3, a conservative es¬ 
timate for high-redshift syste ms consistent with other 
works (Fumagalli et al. 2011a). We also only consider 
the smallest allowed value of Ah i , since this corresponds 
to the most conservative upper limit. The upper limit is 
the highest metallicity for which all model column den¬ 
sities (at the smallest TVhi allowed) were less than the 
measured limits at log U = —3. We refer to these as 
“Type 1” upper limits. Figure [4] is an example of a Type 


1 upper limit. The strongest column density constraint 
for this example came from Si II 1260 due to its large 
oscillator strength, although the Lya forest made it un¬ 
available for many of our systems. 

For absorbers with some detections but not enough 
to fully constrain the posterior distribution in [M/H]- 
log£7 space, the metallicity upper limit is derived from 
the posterior distribution. We take the upper limit to be 
95th percentile of the posterior metallicity distribution. 
We refer to these as “Type 2” upper limits. J080853 is 
presented as an example in the top right of Figure [6j 
Note that the posterior distribution for log U does not 
extend below —3—if it did, then we would reclassify this 
system as a Type 1 upper limit. 

Type 2 limits generally result in lower metallicity lim¬ 
its than Type 1, since measured metal column densities 
are able to constrain the ionization parameter to a larger 
limiting value. Since the posterior distribution for these 
systems includes the largest ionization parameters and 
lowest metallicities, the exact value taken as the limit de¬ 
pends on the range of [M/H] allowed by our priors, but 
in practice our prior distribution was realistic enough 
that large modifications are not physically motivated, 
and small changes to the priors have negligible effect on 
the result. 


3.3.4. Applicability of Single Cloud Models 

Our analysis treated each absorption system as a sin¬ 
gle cloud comprised of gas with minimal phase struc¬ 
ture, which is insufficient t o capture the fu l l stru cture 
and conditions of the CGM. Churchill et al. (2015) con¬ 
structed mock absorption lines through the CGM of a 
simulated z = 0.54 dwarf galaxy, investigating the kine¬ 
matics and phase structure of the gas. For low-ionization 
metals (e.g., Mg II AA2796,2803) as well as H I, ab¬ 
sorption along their simulated sightlines was generally 
dominated by a narrow, single-phase cloud that could be 
readily modeled with typical ionization correction tech¬ 
niques. However, absorption from high-ionization gas 
(e.g., O VI AA1031,1036, C IV) often came from more 
extended structures with varying gas properties, as well 
as gas unassociated with the H I but with a coincident 
line-of-sight velocity. 

This suggests that the influence of C IV and Si IV on 
our results needs to be examined, to gauge whether or 
not they adversely affect our analysis. If we assume that 
an appreciable fraction of the high-ionization metal ab¬ 
sorption along a sightline were due to gas more extended 
than or disjoint from the H I, then the measured high- 
ionization metal column densities would be upper limits 
to the associated H I column densities. 

To assess how this would influence our ionization mod¬ 
eling, we use the LLS along the sightline to J123525 (Fig¬ 
ure pi) as an example. If we treat the C IV and Si IV as 
upper limits, then the only ion with a measured column 
density is C II. Performing the ionization modeling un¬ 
der this assumption we found a posterior distribution 
essentially inverse that of a Type 2 upper limit (Figure 
[6j bottom right): the metallicity is constrained, but the 
ionization parameter is not. This could be treated as a 
Type 1 upper limit, which would have [M/H] < —1.5, 
much larger than the metallicity measured treating C IV 
and Si IV as detections. The explanation for the large 
change can be seen in Figure [5| As discussed in Section 
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Fig. 6.— Top Left Posterior distribution for ionization modeling of J123525. Histograms of [M/H], log U, and Nmare along the 
diagonal, with dashed lines showing the 16th, 50th, and 84th percentiles. The other three plots show the positions of the walkers at each 
step, outlining correlations between the three variables. [M/H] and log U are not strongly dependent on TVhi, but [M/H] and log U are 
tightly correlated. The width in the [M/H]— Nm profile reflects the range of acceptable ionization parameters. The bias in jVhi results 
from moderately larger likelihoods for the metal column densities at larger N-r i values for this system, but does not appreciably influence 
the results. Top Right Same as top left, for J080853. This system gives a Type 2 upper limit: there are upper limit constraints on the 
metallicity (i.e., lower limit constraints on the ionization parameter), but the walkers do not converge to a solution. The dashed lines 
correspond to the 90th, 95th, and 99th percentiles. The red lines are the 95th percentile in [M/H] and the 5th percentile in log U (since 
U and [M/H] are inversely related). Bottom Left J123525 (same as top left), but only half the C IV and Si IV column density, to test 
the effect of interloping gas that is coincident but unaffiliated with the LLS. The posterior is largely similar to that using the measured 
column densities. Bottom Right J123525 again, but with C IV and Si IV treate d as upper limits. In this scenario, the analysis yields a 
metallicity upper limit. Figures formatted with the Python module Triangle ( |Foreman-Mackey et al .|2014 ) 


3.3. H model C II column density is moderately depen¬ 
dent on /Vhi, so without any other lines constraining the 
result, taking the metallicity limit at the smallest viable 
value of A^h i=17.8 and log U =-3 gives a high-metallicity 
upper limit. When all metal lines are treated as detec¬ 
tions, the metallicity is much lower at the same TVhi 
because other detections limit the possible solutions. 

However, since there is a C II detection, treating this 
as a Type 1 upper limit ignores critical information. An 
actual Type 1 upper limit posterior distribution would 
have both [M/H] and log U extending to [M/H] = —4, 


the low metallicity cutoff of our priors. The C II detec¬ 
tion constrains the metallicities for a given value of log U 
and restricts the maximum value of log U. The poste¬ 
rior distribution for this example gives [M/H] = —2.06, 
fairly close to nominal metallicity found for this system 
in spite of the low values of the ionization parameter that 
are included in the posterior. Restricting the posterior 
to logC7 >-3, the modeling gives [M/H]=-2.17. Hence, 
treating C IV and Si IV as upper limits for this LLS 
change the posterior distribution for log J7, but [M/H] 
only changes by several tenths of a dex. 
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Fig. 7.— Metallicities with a bin size of 0.25 dex centered on half¬ 
integer values. The shaded region consists of detections, whereas 
the unshaded portions correspond to upper limits, except for lower 
limit in the metal-blind sample above [M/H]=-l. Although the 
distribution of the detections is not markedly different, the metal- 
poor sample contains far more metallicity upper limits. 


Fig. 8. — Metallicity as a function of redshift for the metal-poor 
(black points) and the metal-blind (green) samples. We maintain 
this color scheme throughout. Arrows indicate upper limits The 
blue point is the IGM, with shaded regions showing the la and 2 a 
limits. Both samples are consistent with having a fraction of their 
metallicities drawn from the IGM. 


Even if some of the high-ionization metals are from in¬ 
terloping gas, ionization modeling still predicts a sizable 
column density from the LLSs. If, for J123525, instead of 
treating the C IV and Si IV as upper limits, we take half 
of the observed column density to be from the LLS, the 
MCMC posterior distribution (Figure [6| bottom left) is 
very similar to that found using the actual measured val¬ 
ues, with a median metallicity of [M/H] = —2.49, slightly 
lower than the metallicity we measure without altering 
the measured column densities. 

The results found for this example generally extend to 
other systems with both high- and low-ionization detec¬ 
tions: treating high ions as detections leads to the inclu¬ 
sion of low log U values in the posterior distribution, but 
ultimately only influences the metallicity by less than 0.5 
dex. 

There are several LLSs for which we detect only high- 
ionization metals, most of which are Type 1 or 2 upper 
limits. Since C II and Si II non-detections tend to place 
strong constrains on metallicities, treating C IV and Si IV 
detections as upper limits for these systems tends to pro¬ 
duce metallicity upper limits within 0.5 dex of the mea¬ 
sured limits. There is one outlier, J115321, a Type 2 up¬ 
per limit that would become a Type 1 upper limit with 
a limiting metallicity that is 1 dex larger. There are two 
LLSs in the metal-blind sample (J234466 & J1025909) 
with metallicity solutions, but only C IV and Si IV de¬ 
tections. We found these systems have posterior distribu¬ 
tions highly constrained by Si II and C II non-detections, 
such that treating the high-ionization detections as up¬ 
per limits produces Type 1 metallicity upper limits that 
are only ^0.5 dex larger than the measured metallicities. 

4. RESULTS 

In Table [2] we list the properties of the quasars and 
LLSs comprising our metal-poor and metal-blind sam¬ 
ples. Of the 17 metal-poor LLSs, nine have metallicity 
upper limits—six Type 1 limits derived from five LLSs 
with no metal detections along with one LLS with only 
C IV, and three Type 2 limits from systems without in¬ 
formation enough information for the MCMC walkers to 


converge. In addition to having generally higher metal- 
licities, the metal-blind survey has only one upper limit 
(Type 2). The metal-blind sample also has one metallic- 
ity lower limit, derived from several ions with saturated 
absorption profiles, with a metallicity well above any seen 
in the metal-poor sample. 

In Figures [7] and [8j we display a histogram of the 
metallicities and a comparison of the LLS metallicities 
with the IGM metallicity as a function of redshift. For 
both the metal-poor and metal-blind samples, measured 
metallicities range from 2.8 to ^—1.8. Only consider¬ 
ing detections, the metal-poor(blind) sample has median 
metallicity of -2.21 (-2.13). Although the distributions 
for metallicities we were able to measure are not strik¬ 
ingly different, the metal-poor sample has a much larger 
fraction of systems for which we were only able to assign 
an upper l imit to the metallicity. 


Simcoe (2011) find at z=2. 4(4.3), the IGM has an 
abundance of [C/H] =-2.90(-3.55) ±0.8 dex, as indicated 
by the shaded region in Figure [5] All of the LLSs in both 
our samples have metallicities within 2a of the diffuse 
IGM metallicity, and several of the systems are in very 
good agreement, having [M/H] ^ —3.0 at z ~ 3.5. Con¬ 
sidering the large fraction of the metal-poor sample con¬ 
stituting metallicity upper limits, this suggests that the 
gas comprising a significant fraction of these absorbers 
has not cycled through a galaxy; if the corresponding 
absorbers represent circumgalactic material, they would 
likely be accreting onto the galactic disk rather than be¬ 
ing expelled. 

We also note that the detections and upper limits are 
not well stratified—although one might expect upper 
limits to fall below the detections, many of the mea¬ 
sured metallicities are below the upper limits. This re¬ 
sults from detections constraining the gas as more highly 
ionized (which typically corresponds to a lower metallic¬ 
ity for the ions we consider) and could mean systems with 
upper limits likely have even smaller metallicities. Addi¬ 
tionally, measuring Type I upper limits at the smallest 
viable value of TVh i and log U results in conservative lim¬ 
its. If we instead use either the median value of allowed 
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TABLE 2 

Lyman Limit Systems 


QSO 


ZQSO 

2 +LS 

An b 

(kms -1 ) 

log iV H 1 
(log[cm —2 ]) 

log U 

n H 

( 10 - 3 cm —3 ) 

t 

(kpc) 

[M/H] d 


Metal-Poor Sample 

J080853—070940 

19.3 

3.841 

3.545 

+80,-80 

17.8-18.6 

> -2.57 

< 3.7 

> 28 

< -2.40 2 

J083832+200142 

18.2 

3.876 

3.476 

+220,-130 

17.9-18.6 

+0.04 
Z.OO — 0 05 

2.8-4.7 

30-40 

— 1 70+ 0 - 09 
i ‘‘ u -0.07 

J085944+212511 

18.8 

3.699 

3.263 

+190,-110 

18.8-19.6 

-2.38 + 0.04 

1.9-3.1 

90-110 

— 1 99 + 0 - 06 
i.yy_Q 07 

J110657+081643 

19.1 

4.268 

4.105 


17.3-17.6 


< 10.9 

> 1.4 

< -2.19 1 

J115321+101112 

19.1 

4.127 

4.038 

+70,-70 

17.7-19.0 

> - 2.22 

< 1.8 

130-190 

< -2.90 2 

J123525+014945 

19.1 

4.031 

3.891 

+150,-150 

17.8-19.0 

- 2.11 + 0.06 

1 . 0 - 1.8 

260-330 

9 07 +O.II 
- 0.10 

J124957—015928 

17.6 

3.638 

3.524 


17.8-19.0 


< 10.1 

> 3.4 

< —2.70 1,d 

J125949+162005 

19.0 

3.707 

3.547 

+ 110,-110 

17.8-18.6 

— 1 Q4 + 0 - 07 
i - b '^-0.09 

0 . 6 - 1.2 

650-890 

-2.92+11 

J130452+023924 

18.4 

3.651 

3.336 

+ 110,-110 

17.9-18.7 

-2.08 + 0.09 

0 . 8 - 1 .7 

360-460 

— 2 .S 1 ++ 

J130452+023924 



3.324 

+90,-90 

17.9-18.6 

> -2.29 

< 2.0 

130-160 

< -3.08 2 

J131056+105530 

19.0 

4.461 

4.200 

+100,+130 

18.2-19.0 

-2.45 + 0.05 

2.6-4.1 

50-60 

9 oq + 0.09 

z -oy _ 0 10 

J134723+002158 

19.3 

4.308 

4.229 

+ 100,-100 

17.5-18.1 

—2 24 +0 ‘° 4 
Z - Z4 -0.05 

1 . 6 - 2 .5 

80-140 

— 2 05 4-0 ' 10 
z - uo — 0.08 

J144027+173038 

19.5 

3.674 

3.566 


19.2-19.6 


< 10.0 

> 5.2 

< - 1 . 68 1 

J144405+165621 

18.9 

3.745 

3.551 

+150,-130 

17.3-17.6 

-2.05 + 0.12 

1 . 1 - 1 .7 

170-290 

—2 4i +°. 20 
-^-0.23 

J144405+165621 



3.471 


17.5-17.8 


< 10.1 

> 2.2 

< - 2 . 21 1 

J155255+145432 

19.9 

4.105 

3.954 


17.8-18.4 


< 10.5 

> 3.3 

< — 2.03 1 

J160320+072104 

19.5 

4.393 

4.375 

+60,-60 

17.8-18.4 


< 11.8 

> 2.8 

< — 1.92 1 

Metal-Blind Sample 

J001022-003701 

18.4 

3.153 

3.116 

- 120,+120 

17.8-18.8 

9 11+0.04 

Z- ^ ^ — 0.05 

1 . 0 - 1 .7 

270-400 

—2 19+ 0 - 07 

J014850-090712 

18.0 

3.322 

2.996 

-110,+170 

17.8-18.7 

—2.46 + 0.03 

2.4-3 .8 

40-70 

— 2.06 + 0.05 

J030341-002321 

17.7 

3.229 

2.941 

-150,+150 

18.7-19.2 

—2.14 + 0.02 

1 . 2 - 1.8 

300-320 

—2.07 + 0.03 

H0449-1325 


3.107 

2.997 

-70,+70 

17.8-18.2 

> -2.55 

< 3.8 

> 30 

< — 2.69 2 

J093153-000051 

18.7 

3.209 

2.927 

-150,+150 

18.4-19.3 

9 01 +0.03 
Z - Oi -0.04 

1 .8-2.8 

120-140 

9 9 P+O.O 6 
z-zu — 0.05 

J102509+045246 

19.2 

3.243 

3.130 

-70,-70 

17.8-18.7 

9 09 +O.I 3 
z - oz - 0.12 

1.3-3.3 

100-140 

9 715 +O.I 3 

Z - /O - 0.12 

J161545+060852 

18.2 

3.062 

2.988 

-110,+160 

17.8-19.5 

9 71 +0.06 
z - 1 -*--0.08 

4.1-7.6 

10-20 

—2 02 +0 ‘ 14 
z - uz —0.11 

J223819-092106 

18.0 

3.278 

3.127 

-160,+170 

17.8-19.0 

> -2.60 

< 4.1 

> 21 

> -0.75 

J233446-090812 

18.0 

3.351 

3.226 

-90,+100 

17.8-18.6 

—2 05 4-0 ' 07 
Z - UO -0.06 

0 . 8 - 1 .5 

380-520 

9 71 + 0.11 
Z -' 1 -0.15 

UM184 


3.021 

2.929 

-110,+110 

18.5-19.2 

9 151 +0.01 
Z - Oi -0.02 

2.9-4.3 

40-50 

1 70 +O.O 3 

4 '‘ ° —0.05 


a Quasar r-band magnitude 

b Velocity width of absorbers, based on metal detections. Dots denote that all lines had no absorption. 
c Superscripts indicate Type 1 and Type 2 limits 

d Using column density limits from the MIKE+ FIRE spectra, we find [M/H] < —2.90. 


TVhi or the median log U from systems with detections, 
Type 1 limits become stricter by up to ^ 0.5 dex. 

4.1. Ionization Parameters 

Larger values of the ionization parameter U gener¬ 
ally correspond to lower derived abundances for the 
ions we considered, so systematically overestimated ion¬ 
ization parameters wou ld drive our res u lts to artifi¬ 
cially low metallicities. Fumagalli et al. (2011a) com¬ 
piled results from published high-redshift LLSs (2 >1.5, 
their Figure S6) and determine that all systems have 
—3 < log U < —1, with most systems in the range 
of —3 < logC7 < —2.5. These systems generally have 
smaller redshifts than our sample. For all of the ab¬ 
sorbers with detected heavy-element absorption in our 
metal-poor sample, we find —2.7 < log U < —1.9, with 
most clustered around log U = — 2.3 (Figure [ 9 ]). 

There are several plausible explanations for the mod¬ 
erately larger ionization parameters we m e asured . The 
ionization parameters in Fumagalli et al. (2011a) were 
compiled from studies using varying modeling techniques 


and ionizing radiation spectra. In addition, the complete 
literature sample comprises a relatively small sample, 
roughly the same size as our survey, often with redshifts 
and H I column densities different than our sample, sug¬ 
gesting that they may not form a unifor mly selecte d com ¬ 
parison group for our results. Indeed, Fuma galli et al. 
(2011a) used the values from the literature only to show 
that log U = —3 is a justifiable lower limit. The ab¬ 
sorbers in Crighton et al. (2015) are at z = 2.5, and have 
sub-LLSs H I column densities. 


Lehner et al. (2013) studied absorbers with 16.2 < 
log 7 Vht < 1 8.5 at 2 ^ 1 and found log U = —3.3 =b 0.6. 


Werk et al.| (|2014|) looked at LLSs located near L * galax¬ 
ies at 2 5: 0.2 (overlapping the sample from Lehner 


et al. 2013) and measured a mean ionization parame 
ter log ~U = —2.8. A trend of LLS ionization parameter 
moderately increasing with redshift is supported by the 
studies previously mentioned and our sample, although 
the total sample size is small. 

The hydrogen density, tir, for each LLS follows from 
the ionization parameter U from Equation^ To calcu- 
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Fig. 9. — Ionization parameters derived for both of our samples 
and from the literature. In our work, black and green correspond 
to the metal-blind and metal-poor samples, respectively. LLSs giv¬ 
ing metallicity limits are not included. In the litera ture histogram, 
the black corresponds to |Fumagalli et al. | (|2011a| and references 
therein ), and the hashed region is absorbers trom JCrighton et al.| 
(2015|), who perform an analysis similar to ours. The literature 
sample is generally at lower redshift,and contains numerous sys¬ 
tems with A^hi lower than our sample. 


late <f>, we used Equation [2] using the sha pe of the UVB 
spectrum from |Haardt & Vladau (2012), renormalized 
to match the observed H I photoionization rates from 
(2013). Under the uniformity assump- 
the gas clouds can be estimated via 


Becker & Bolton 
tion, the sizes o 


£ = A^hi/(xhi^h) ? where the neutral hydrogen ioniza¬ 
tion fraction %hi is output by the ionization models. In 
Table [2] we list the range of densities corresponding to 
the ionization parameter posterior distribution, and the 
range of sizes that follow assuming the central value of 
the log U over the range of viable TVhi- Typical densities 
are of order 10 -3 cm -3 , and typical lengths are a few tens 
to hundreds of kiloparsecs, with a median of 160 kpc. 

In Figure 10, we plot [M/H] over the overdensities 
S = nn/nH of the absorbers, where % is the cosmic 
mean baryon density at the redshift of the absorber: 
n h = UfrPcfl + z) 3 /m p , where 1^5=0.04 is the cosmic 
baryon density relative to p c , the critical density of the 
universe, and m p is the proton mass. The systems 
with detections have typical overdensities ranging from 
5 ~ 10-100. This figure also portrays how conservative 
the limits are: the Type 1 upper limits correspond to 
limiting overdensities of at least two-fold greater than 
the detections, implying the metallicities may be signifi¬ 
cantly lower than the limits we adopted. 


4.2. Effect of High-Resolution and Infrared Spectrum 

The high-resolution MIKE spectrum of J124957 con¬ 
firmed the lack of metal lines in the MagE spectrum of 
the same object was not a resolution effect, and the IR 
FIRE spectrum likewise shows no absorption at expected 
locations. In Figure [TT] we compare the SDSS, MagE, and 
MIKE spectra for J124957, as well as showing cutouts of 
several portions of the FIRE spectrum where absorption 
is expected. It is evident that the four-fold increase in 
resolution provided by MIKE does not reveal any weak 
lines in this particular case. 

We obtained stricter column density upper limits with 
the high-resolution spectrum, allowing us to measure a 



6 

Fig. 10. — Metallicities and overdensities derived from ioniza¬ 
tion modeling. Uncertainties on the ionization parame ter (w hich 
corresp on ds t o nn) and H I photoionization rates from (Becker &| 
Bolton] (2013|) used to renormalize the ionizing spectrum contribute 
roughly equally to the overdensity uncertainties. The large limit¬ 
ing overdensities on the Type 1 upper limits suggest the actual 
metallicities may be much lower. 


metallicity upper limit of [M/H] < —2.90, 0.2 dex less 
than the upper limit measured with MagE. The limits 
we measure from the FIRE spectrum are not as strict, 
since it has lower resolution. When we use column den¬ 
sity limits from both MIKE and FIRE, the ions mea¬ 
sured with FIRE do not influence the result. Using only 
ions from the FIRE spectrum, we find an upper limit of 
[M/H] < -2.60. 

It remains to be seen if high-resolution observations 
could reveal weaker lines in other examples, but in 
light of current sensitivities and to maintain consistency 
among the sample, we use the metallicity limit obtained 
from the MagE spectrum in subsequent analysis. 


4.3. Survival Analysis 

A complete description of the distribution of metallic¬ 
ities needs to incorporate information provided by both 
detections and upper limits. To that end, we employ sur¬ 
vival analysis methods developed to deal with censored 
data sets containing a mixture of measurements and lim¬ 
its. 

For univariate data, the Kaplan-Meier estimator pro¬ 
vides a general, non-parametric maximum-likelihood es¬ 
timate of the population from which a censored sample 
was drawn. For details on the applic atio n of the K aplan- 
Meier method to similar datasets, see Simcoe et al. (2004) 
and references therein. 

Briefly, the Kaplan-Meier method constructs a cumu¬ 
lative distribution function for the sample, handling am¬ 
biguities introduced by upper limits by only including 
them in probability calculations when they are guar¬ 
anteed to be unambiguous. For example, an upper 
limit of [M/H] < —2.5 is guaranteed to be less than 
[M/H] = —2.4, but may not be less than [M/H] = —2.6 
and will not be treated as such. The resulting cumu¬ 
lative distribution is a piece-wise function that remains 
constant at upper limits and jumps at detections. 

The Kaplan-Meier method requires the sample to sat¬ 
isfy two criteria. First, the upper limits must be in¬ 
dependent. This is clearly true for our sample, where 
each measurement is drawn from a different absorber and 
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Fig. 11. — Top: Comparison of SDSS, MagE, and MIKE data over a portion of the J124957 spectrum. The MIKE spectrum has a blue 
line indicating the transition from the blue and red arms of the instrument. Bottom: Portions of the normalized spectra around expected 
LLS lines, showing no absorption. The absorption offset by 250kms —1 in the Si IV panel is an interloping C IV line at z = 3.102. To the 
right are three cutouts of the FIRE spectrum, also showing no absorption. 


most are from different sightlines. Second, the probabil¬ 
ity of a measurement being censored must be uncorre¬ 
lated with the value of the measurement itself. While 
our sample likely does not strictly meet this criterion, 
since lower metallicity systems are more likely to result 
in non-detections, there is a characteristic to our sur¬ 
vey that preserves the randomness of the censoring: all 
targets were selected using the same criterion, so the pri¬ 
ors on metallicity are uniform across the sample. Hence, 
the selection method should adequately randomize the 
censoring. Also, since an LLS observation resulting in a 
metallicity limit partially depends on the signal-to-noise 
ratio of the spectrum, the brightness of the quasar and 
observing conditions also serv e as r andomizing factors. 
From the discussion in Section |3.3.1[ metal column den¬ 
sities generally do not depend strongly on TVhi, so TVhi 
should not bias whether or not a system gives a metallic- 
ity upper limit. Since the metallicity measurements and 
limits we find are not segregated such that all limits fall 
below detections, the Kaplan-Meier method is applica¬ 
ble. 

We calculate the Kaplan-Meie r distributions for our 


samples using ASUR Y Rev 1.2 (Isobe & Feigelson 


Lavalley et al 1 19921), which implements the methods dis- 


1990 


cussed in Feigelson Kelson (1985). The resulting cu¬ 
mulative distributions are shown in Figure [12] 

We also extrapolate our results for eachsample in¬ 
dependently to estimate the entire LLS population of 
z ~ 3.73, the mean redshift of the metal-poor sample. 
Since the metal-blind sample is at lower redshift and 
metallicities evolve with redshift, we can apply a shift 


to the entire cumulative distribution function (CDF) to 
bring it to the same redshift as the metal-poor sample. 
Taking the mean slope of the IGM and D LA m etallicity 
with redshift, [M/H] oc 0.2 8z (see Section 5.1), we shift 
the metal-blind CDF by the difference between the mean 
redshifts of the samples. 

Alternatively, we may estimate the full LLS CDF from 
the metal-poor sample if we assume all systems with 
metal lines detected in SDSS have higher metallicity than 
those that do not. This assertion is demonstrably false 
for some individual cases, but lacking the (forthcoming) 
fully unbiased set of H I-selected LLS metallicities, it can 
represent a first attempt at generalization. In this case 
we may construct the CDF for the entire SDSS sample 
as: 


PsDSS = xP-MagE + (1 ~ x) 

Number of SDSS LLSs meeting metal-blind criteria 
1 ~ Number of SDSS LLSs 

In other words, we scale the CDF by the fraction of 
LLSs matching our selection criterion (x), then add to it 
the fraction of LLSs that do not. Since, under our as¬ 
sumption, all LLSs not meeting our criteria have metal¬ 
licities larger than those that do, this approximates the 
CDF for metallicities in the range probed by our metal- 
poor sample. In our scan of the SDSS spectra, we found 
x = 0.48. We stress that variations in signal-to-noise and 
ionization fractions may result in some SDSS LLSs that 
failed to meet our initial selection criteria having lower 
metallicity than others that did. 

In Figure [12] we compare these extrapolations to the 
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Fig. 12. — Left: Cumulative distribution functions for the 
metal-poor (black) and metal-blind (green) LLS samples, con¬ 
structed using the Kaplan-Meier Estimator. Right: Estimates 
of the full LLS CDF at z ~ 3.7 extrapolated from the metal-poor 
(black-dashed) and metal-blind (green-dashed) samples. They dis¬ 
agree at high metallicity since the metal-poor sample does not 
probe this region. 


metal-poor CDF and to each other. They diverge at high 
metallicity, because the assumptions on the extrapola¬ 
tion from the metal-poor sample place an unrealistic floor 
on the corresponding CDF, which is most relevant at 
higher metallicity. They are in fairly good agreement for 
metallicities below [M/H] = —2.5; both extrapolations 
suggest ^ 20% of LLS at 2 ~ 3.73 have [M/H] < —2.5, a 
value roughly la above the measured IGM abundance. 

5. DISCUSSION 

In this section, we leverage our dataset to extract phys¬ 
ical and cosmological details concerning low-metallicity 
LLSs and compare with expectations for cold flows from 
simulations. Since our sample size is small, our inten¬ 
tion is to establish an order-of-magnitude for several key 
properties. 


5.1. Interpreting the Distribution in the Context of 
Cold-Flows 

In order to assess whether or not low-metallicity LLSs 
are consistent with being cold flows, we compare the 
metallicity CDFs of both the metal-poor and metal-blind 
samples to a toy model parent CDF. Motivated by the 
bimo dal me tallicity distribution found at low redshift 
(Lehner et al.|2013), we assume a mixed Gaussian model 
where the absorbers are drawn from a combination of two 
different parent populations, one being the IGM (repre¬ 
senting potential accretion flows) and the other having 
more highly enriched gas that has been polluted by a lo¬ 
cal host galaxy. We refer to absorbers drawn from the 
IGM distribution as cold-flow candidates (CFCs). 

We assume the parent distribution in [M/H] for CFCs 
is the same as the IGM’s, which we interpolate from the 
measurements of Simcoe (2011) to the mean redshift of 
the sample. Note that this study used the same ionizing 
background spectrum to measure metallicities, so any 
systematics from uncertainty in a particular realization 
of the UVB spectrum are common to both studies. 

The parent [M/H] distribution of enriched CGM gas is 
less clear, but for this study we associate this phase with 
DLA abundances. The exact physical structures giving 
rise to DLAs are neither fully understood nor expected to 
be uniform, but they are thought to be locally enriched. 
Since DLAs have systematically lower abund ances th an 
H II regions measured in emission (e.g., Sanders et al 


2015), this is a conservative choice to represent the non- 
CFC branch (using larger metallicities for this branch 
would give a larger fraction of CFCs). We model enriched 
gas using a lognormal distribution with parameters given 
by DLA measurements, since DL As are thought to orig - 


inate from gas in galaxies (e.g., Rafelski et al. 2011) 


We use DLA metallicities from^iafelski et al. (|2012| and 
references therein). Compared to LLSs, DLA metallici¬ 
ties have been analyzed more extensively since DLAs are 
predominantly neutral and tend to have small ionization 
corrections that do not require modeling. In addition, 
the Ly<a damping wings allow for accurate H I column 
densities in moderate-resolution spectra. 

To compare with the metal-poor (blind) sample we av¬ 
erage over all DLAs between z = 3.26-4.37 (2.90-3.25). 
The model metallicity probability distribution is then 


p([ M/H]) = /igm Pigm([M/H]) + (1 —/igm)pdla([M/H]) 

where Pigm and Pdla are Gaussian metallicity distribu¬ 
tions with parameters summarized in Table i and /igm 
is the fraction of model distribution drawn from the IGM 
that we are estimating (i.e., the fraction of CFCs). 

We performed Kolmogorov-Smirnov tests to find which 
IGM fractions are allowed by the two measured distribu¬ 
tions (Figure [l3| ). We list the values of /igm allowed 
within 68 and 95% confidence in Table [3] 

Performing a least-squares fit, we find a best fit to the 
CDFs with /igm = 0.71 and 0.48 for the metal-poor and 
metal-blind samples, respectively. We tested several dif¬ 
ferent likelihood functions and found approximately the 
same maximum likelihood values for /igm- We adopt 
these best-fit /igm values, with errors given by the 68% 
confidence intervals for the remainder of the paper. We 
caution that the small sample sizes enable intrinsic vari¬ 
ation within the LLS population to appreciably influence 
the results. 

Since 48% of the SDSS sample meets our metal-poor 
criteria, assuming only systems passing our initial cuts 
can be cold-flow candidates implies that the range of ac¬ 
ceptable values for /igm for the entire z ^ 3.7 LLS pop¬ 
ulation is 0.34 ± 0.06. This is somewhat less than /igm 
for the metal-blind sample, suggesting that the assump¬ 
tion regarding the S DSS sample may be questionable, as 
discussed in Section [T3j although sample variance may 
also account for the disagreement. 

In Figure [l4j we show the best-fit CDF and 68% con¬ 
fidence intervals for both samples, as well as the best-fit 
probability distribution function with the relative contri¬ 
butions from the enriched and unenriched parent popu¬ 
lations. 


5.2. Mass Fraction of Candidate Cold Flows 

Of particular interest is Ucfc : the mass fraction of the 
Universe (relative to the critical density) at z ~ 4 con¬ 
tained in LLSs with IGM metallicities (i.e., CFCs). This 
quantity can be compared to simulations to test whether 
the global mass contained in cold flows agrees. Using our 
measurements and ionization models we make an order- 
of-magnitude estimate of this quantity for comparison 
with simulations. 

The ratio of the CFC mass density to the cosmological 
critical density, p c , is given by: 

Dcfc = ——— [ NH(Nm)fcFc{Niu)dNm 
C Pc Jo 
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TABLE 3 

Metallicity Distribution Parameters 



z 

[M/H]igm 

<TlGM 

[M/H] DLA 

^DLA 

/igm 

68% c.i. 

95% c.i. 

Metal-poor 

3.73 

-3.36 

0.8 

-1.69 

0.48 

0.71 

0.60-0.84 

0.44-0.96 

Metal-blind 

3.04 

-3.12 

0.8 

-1.49 

0.52 

0.48 

0.34-0.62 

0.18-0.79 



f\GM 


Fig. 13 . — Kolmogorov-Smirnov test P-values as a function of 
flGM.foY the metal-poor (black) and metal-blind (green) samples. 
The null hypothesis is that the observed CDF is drawn from a 
parent sample having a fraction /igm of its metallicities coming 
from the IGM distribution. Thick bars denote the 68 and 95% 
confidence intervals 


where Ah (Am) is the total hydrogen column density of 
an absorber, /cfc(^hi) is the frequency distribution of 
CFCs as a function of neutral hydrogen column density, 
/r is the reduced mass of the gas and ran is the mass of 
the hydrogen atom. 

To computethe integral, we need to assume a form for 
/cfc(Nhi) POW10 find /(A H i) for 2? = 3.7 LLSs can be 
fit by: 


/lls(A H i) = 


IO-^ss^-o.s 17,5 < l 0 gA^ HI < 19 
lO 2 - 75 ^! 1 - 2 19 < log Vhi < 20.3 


To estimate /cFc(Nfin), we multiply /lls by the frac¬ 
tion of LLSs that are cold-flow candidates, /igm- (Note 
/igm is the fraction of LLS that are CFCs and /cfc 
is the frequency distribution of th ese LLSs). We take 


/igm=0.34±0.06 from Section 5.1 


We assume a reduced mass fi = 1.3, appropriate for 
absorbers with 75% H and 25% He by mass. The to¬ 
tal hydrogen column density is readily found from the 
H I column density Ah = (^hA&hi)-Nhi, with the ratio 
of ionized-to-total hydrogen coming from the ionization 
model. We use the median value derived from ionization 
solutions in both of our samples: (nn/^Hi) = 0.0063. 
The median value for each sample separately is similar. 
Due to our small sample size, we cannot adequately mea¬ 
sure and include in the calculation any variation in /cfc 
and (nn/^Hi) with H I column density. 

We also need to set bounds on the integral. For the 
lower bound, we use log Am = 17.5, t he stated sensitiv¬ 
ity limit of the LLS survey in PQW10[ Although lower 
column density absorbers are more numerous, larger col¬ 
umn density absorbers dominate the mass density, so 
the result is largely insensitive to the lower bound. As 



Fig. 14. — Left: Model CDFs corresponding to the best-fit /igm 
(blue) overlaid onto measured CDFs. The shaded region encom¬ 
passes 68 % confidence on /igm- Right: PDF corresponding to 
the best-fit /igm and the contributions from the IGM (dotted) 
and DLA (dashed) distributions. 


an upper bound we take log Am = 19.5, a typical value 
where systems transition from being considered LLSs to 
sub-DLAs. 

With these parameters, we obtain Hcfc = 0.0017. 
Comparing this to the cosmic baryon density, fib = 0.04, 
we find roughly 5% of baryons at z ~ 3.7 are contained 
in cold-flow candidate LLSs. Note this calculation is sen¬ 
sitive to both the maximum column density used in in¬ 
tegration and the H I ionization fraction; increasing (de¬ 
creasing) the maximum Am by 0.2 dex increases (de¬ 
creases) the result by a factor of ~ 1.5, and the H I 
ionization fractions for individual systems can vary from 
the median by a factor of ~5. As we discuss in Section 
|5.4[ our result is fairly consistent with simulated results. 


5.3. Comparison with Other Observations 

In Figure [l5j we compare th e metallicities of our sam- 
ples to DLA metallicities from Rafelski et al. (2012, and 
references therein). DLAs generally have higher metallic¬ 
ities than both our metal-poor and metal-blind samples, 
and there are suggestions that DLAs have a metallicity 
floor, which many LLSs are below. It is clear from this 
comparison that LLSs and DLAs at high redshift differ 
significantly in their metallicity distributions._ 

Also shown are low-redshift LLSs from ILehner et al.l 
(2013, and references therein). They categorized their 
H I systems as LLS (16.2 < log Am < 19) or super 
LLS (SLLS, 19 < log Ah i < 20.3). Their LLS sample 
is mostly composed of systems with log Ah i < 17.5, the 
cutoff for our survey, so some differences are expected. 

The low-redshift LLS population clearly exhibits a 
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Fig. 15 . — Comparison of our results with DLAs (orange) and 
LLSs (various) from the literature. At high redshift, the DLA 
and LLS metallicity distributions are more clearly different than at 
low redshift. The bimodality seen in low-redshift LLS metallicities 
is not evident at high redshift, although the population of high- 
redshift LLSs has not been fully explored. 


metallicity bimodality, with most of the lower-metallicity 
branch below most DLAs, although the difference be¬ 
tween the LLS and DLA populations is not as empha¬ 
sized as at higher redshift. While a bimodal model fits 
our high-redshift sample well, the two populations blend 
together more sm oothly than they do at lower redshift, 
where [Lehner et al.| see very few systems at intermedi¬ 
ate abundances. This may be a result of many of the 
lower abundances at high-redshift producing upper lim¬ 
its rather than measurements. 

Several LLSs drawn from the literature (Table [I]) that 
exhibit low metallicities and are claimed as potential ev¬ 
idence of cold flows are also included in Figure [15] Our 
work corroborates the finding of low-metallicity LLSs and 
provides some statistical context for the population from 
which they are draw n. The two high-redshi ft metallicity 
upper limits are from|Fumagalli et al.|(|2013|). Using high- 
resolution spectra, they were able to model and subtract 
the Ly<a forest to obtain column-density upper limits for 
C III A977 and Si III A1206, which provide tighter con¬ 
straints than the ions available in our medium-resolution 
spectra (see their Figure S5). 


5.4. Comparison with Simulations 

By comparing our sample with structures having anal¬ 
ogous properties in simulations, we explore the agree¬ 
ment between simulations and observations_and gain in- 
sight into the nature of metal-poor LLSs. Fumagalli et al. 
(|2011b) simulated absorption profiles produced by cold 
hows at z ~ 2.3 and found that much of the gas is ion¬ 
ized by the UV background, appearing mostly as LLSs. 
They determined that DLAs have higher metallicities 
than LLSs and SLLSs, with DLA metallicities fairly con¬ 
sistent with observed systems, and the authors suggested 
metal-poor LLSs may therefore be an observational sig¬ 
nature of cold-how accretion. Their simulations predict 
that most of the cold-how observational signatures are 
LLSs with 17 < logTVni < 18, with a peak metallicity 


TABLE 4 

Low-Metallicity LLSs in Literature 


QSO 

2 +LS 

[M/H] 

Source 



PG1630+377 

0.274 

-1.71 + 0.06 | 

|Lehner et al. (2013 


J144535+291905 

2.44 

-2.0 + 0.17 

CrigRton et al. (2UJ 

L3 


HE 0940-1050 

2.917 

-2.93 + 0.13 

lEevsUakov et al. (2UU 

3) 

J113418+574204 

3.411 

< -4.2 

Fumagalli et al. ( 

2013 


Q0956+122 

3.096 

< -3.8 

Fumagalli et al. ( 

2013 



for cold stream LLSs of one-hundredth solar. 

While our measured metallicities tend to be somewhat 
lower, enrichment of the IGM between z = 3.7 and 
z = 2.3 may account for some of the difference. An 
additional difference is that our metal-poor sample is 
mostly composed of absorbers with 18 < logTVni < 19, 
the higher end of the LLS column-density distrib ution. 

A l i mitatio n of the simulations employed in |Fuma galli 
et al. ]2011 b| and s i milar s tudies is the simulated volume. 


In 


alii et al. (2011b), zoom-in simulations of seven 
halos were considered. This restricts analysis of the cov¬ 
ering fraction of LLSs to roughly a few times the virial 
radius of each galaxy. Although these simulations al¬ 
lowed for a descriptive picture of the neutral-gas content 
immediately around the seven relatively massive galaxies 
presented, it remains unclear if a random quasar sight line 
is more likely to intersect a LLS in the immediate vicinity 
of one of these systems or elsewhere in the cosmic web. 
Such questions can only be addressed with larger sim¬ 
ulation volumes. F ull volume cosmological simulations 
can complement the Fumagalli et al. (2011b) analysis by 
probing the full range of LLSs that are probed through 
quasar-selected samples. 

To demonstrate this point, we briefly consider the 
global distribution of neutral hydroge n in t h e full- volume 
cosmological simulations presented in Bird et al. (2014). 
These simulations were run using the cosmological hy- 
drodynamical simulation code AREP0 (Springel 2010) in 
a periodic box of size 10 h -1 Mpc. The simulations con¬ 
tain 512 3 dark matter particles and a similar number of 
baryon resolution elements yielding a mass resolution of 
1.4 x 10 5 Mq —ab out an order-of-magnitu de larger than 
that presented in Fumagalli et al.| (201 lb|) . The simula- 
tion physics is the same as in in the lllustris Simulation 
(jVogelsberger et al. 2014) which importantly includes 
star-formation driven winds at a level that allows for ap¬ 
propriate evolution of the galaxy stel lar-mass function 
(Torrey et al. 2014 Genel et al. 2014). Neutral hydro¬ 
gen fractions are obtained assuming a u niform UV back- 


ground, with self-shielding corrections (Rahmati et al. 


2013). Complet e simulation and p ost processing details 


are presented in Bird et al. (2014). 

We examine the simulati ons by applying a similar tech¬ 
nique to that used in both Fumagalli et al. (2011b) and 
Bird et al. (2014). Specifically, we project the neutral 


hydrogen and mass-weighted average metallicity onto a 
two-dimensional grid. We do this only along one pro¬ 
jected direction but do not expect our results to be mod¬ 
ified if we considered other projections. We employ a 
grid of 16,000 by 16,000 cells, which results in converged 
neutral-hydrogen column-density distribution functions. 
Projecting the full 10 Mpc box onto a single grid could 
boost the LLS number count by adding multiple, well- 
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Fig. 16. — H I column density (left) and metallicity (right) for sightlines with logA^ni > 17.5 in a 1 Mpc thick slice of our cosmological 
simulation. Lower column density material appears as semi-transparent blue. The simulated LLS gas within halos tends to be enriched, 
with metal-poor LLSs tracing the cosmic web. This suggests either inflowing material is being artificially contaminated by feedback or 
observed metal-poor LLSs are outside of large galactic halos. 


separated low column density systems. To minimize this 
effect, we use ten slices each with a thickness of 1 Mpc. 
We treat each pixel as an independent line of sight. 

In Figure [l6j we show a map of the neutral-hydrogen 
column density through one slice of the simulated box 
at z = 3.5, truncated at logTVni = 17.5 (to empha¬ 
size LLSs), and the corresponding map of metallici- 
ties. Lower column density material appears as semi¬ 
transparent blue in the column-density map. There is 
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Fig. 17. — Cumulative fraction of LLS abundances from cos¬ 
mological simulation at z = 3.5. The metal-blind CDF (shifted to 
z = 3.5 according to the redshift evolution in IGM and DLA metal- 
licities) is in fairly good agreement, while the metal-poor CDF is 
lower, as expected since it is a metallicity biased sample. The 
full LLS distributio n ext rapolated from the metal-poor sample as 
discussed in Section |4.3| is in agreement at low metallicities. 


LLS-level column density material tracing the cosmic 
web extending beyond the virial radii of galaxies in the 
simulated volume. 

Since we do not have explicit projected offsets for the 
observed quasar sight lines, it is hard to know which part 
of the IGM/CGM is being probed. It is possible some 
fraction of our observed CFCs intersect material still in 
the IGM, well outside of the halo virial radius. In those 
cases, it is not immediately clear over what timescale the 
observed neutral gas will fall through the virial radius 
nor whether it will stay neutral (or be shock heated) as 
it is accreted. Additionally, most of the simulated LLSs 
within a halo virial radius are enriched, with the metal- 
poor LLSs tending to trace the cosmic web. This suggests 
that either the feedback prescription in the simulation 
artificially contaminates inflowing material or much of 
our observed sample is inter-halo material, as opposed 
to intra-halo. 

We can directly address global statistics of the simu¬ 
lated LLS population. Adopting standard column den¬ 
sity limits (17.5 < logTVni < 19.5), we find roughly 7% 
of the simulation baryons reside in LLSs, and it has been 
previously shown that the column density distribution 
function for neutral hydrogen in t hese mo d els is r eason- 
ably consistent with observations (Bird et al. 20111 The 
simulated LLS CDF at redshift z = 3.5 is shown in Fig¬ 
ure 17 Our metal-blind sample CDF (when shifted to 
z = 3.5f] is in fairly good agreement with the simu¬ 
lated LLSXDF with a peak metallicity for LLSs around 
one-hundredth solar. Although we are sampling a sig¬ 
nificantly la rger volume, thi s result is similar to that 
presented in |Fumagalli et al.| (|2011b|) . An extrapolation 
of the full LLS population from the metal-blind sam- 


9 The agreement between the unadjusted metal-blind sample 
and the simulation at z = 3 is similar. 
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pie is in agreement at low metallicities, but diverges at 
larger metallicities where the simple extrapolation is in¬ 
adequate. 

We consider all LLSs with [M/H] < —2.5 to be 
CFCs. The fraction of LLSs constituting CFCs is ap¬ 
proximately the ratio of the LLS and CFC mass densi¬ 
ties, flcFc/^LLS = 0.312. This is in accordance with the 
observational result for /igm extrapolated to the entire 
LLS population. 

However, we find the derived Ucfc to be smaller for 
the simulations by a factor of ~ 2. Given that the sim¬ 
ulations Jhave both a similar TVhi distribution function 
(Bird et al. 2014) and metal-poor fraction when com¬ 
pared with observations, it is likely this offset is driven 
by the applied ionization corrections. In our IIcfc cal¬ 
culation based on observations, we assumed the hydro¬ 
gen neutral fraction, nni/nn, is independent of TVhi due 
to an insufficient amount of data to measure any such 
relationship. Since the integrals are heavily weighted to¬ 
wards larger TVhi systems, unaccounted for variations in 
neutral fraction could systematically and significantly in¬ 
fluence the calculation. Hence, the disagreement between 
simulated and observed values of Hcfc and ^lls may be 
resolved by improved observational statistics and does 
not necessarily undermine the agreement on the fraction 
of LLSs that are metal poor. 

Additionally, we suspect the prevalence of metal-poor 
LLSs in the simulations can be influenced by (i) the 
specifics of the adopted feedback model and (ii) the mass 
and spatial resolution. Further investigations on both of 
these fronts are warranted, alongside developing an un¬ 
derstanding of the dependence of the hydrogen neutral 
fraction on TVhi- Full-volume cosmological simulations 
that are able to simultaneously reproduce the TVhi dis¬ 
tribution function as well as the low-metallicity tail of the 
LLS distribution presented in this paper will be helpful in 
identifying the fraction of low-metallicity LLSs residing 
within halo virial radii and, thus, the true mass density 
of cold flows. 


6. SUMMARY 

We have completed a medium-resolution spectroscopic 
survey of 17 Lyman-limit systems exhibiting no statisti¬ 
cally significant metal-line absorption in the SDSS dis¬ 
covery spectra (FWHM « 150kms _1 ) to probe the low 
abundance end of the LLS population. The main results 
are as follows: 

1. Five of the LLSs exhibit no statistically significant 
absorption at any of the available metal transitions 
at MagE resolution (FWHM = 60.7kms _1 ) In to¬ 
tal we found nine metallicity upper limits, rang¬ 
ing from [M/H] < —1.68 to < —3.08, with three 
of the upper limits below [M/H] = —2.50. The 
eight remaining LLSs have metallicities ranging 
from [M/H] = —2 to —3 and ionization parameters 
ranging from logU = —1.9 to —2.6. The median 
metallicity for the detections is [M/H] = —2.21. 

2. A sample of ten LLSs at z « 3 selected blindly with 


respect to metal line absorption exhibits somewhat 
different properties. Although the median for the 
systems with measured metallicities is roughly the 
same at [M/H] = —2.13, only one of the sys¬ 
tems has no metal absorption lines. Addition¬ 
ally, this sample contains one LLS with saturated 
metal lines, leading to a metallicity lower-limit of 
[M/H] > —0.75. Taking into the account that over 
half of the metal-poor sample LLSs have metallic- 
ity upper limits, the two samples may have very 
different metallicity distributions, as demonstrated 
using survival statistics. 

3. LLSs in both samples have typical densities of (1 — 
5) x 10 -3 cm -3 , with corresponding overdensities 
ranging from 20-200. Length scales span several 
tens to hundreds of kiloparsecs, with a median of 
160 kpc. 

4. From the cumulative distribution function that re¬ 
sults from a survival analysis of the detections and 
limits, the metal-poor sample is consistent with 
0.7llo.ii °f the metallicities being drawn from the 
IGM metallicity distribution. Nearly half of the 
LLSs in SDSS spectra meet our criterion for be¬ 
ing metal poor, implying that 28-40% of LLSs 
at z — 3.2-4.4 have IGM-consistent metallicities. 
The metal-blind sample is consistent with an IGM 
metallicity fraction of 0.48 ^q'i 2- A comparison be¬ 
tween LLSs and DLAs shows that have distinct 
metallicity distributions, with many LLSs having 
metallicities below the DLA metallicity floor. 

5. We find the cosmic density of low-metallicity LLSs 
(cold-flow candidates) to be Dcfc ~ 0.0017, ac¬ 
counting for ^ 5% of the total baryonic mass bud¬ 
get at this redshift. This is roughly twice the 
baryonic fraction of CFCs in simulations, with the 
disagreement likely attributable to limited infor¬ 
mation of the hydrogen neutral fraction and fre¬ 
quency distribution of low-metallicity LLSs. Sim¬ 
ulations agree with our observed fraction of metal- 
poor LLSs, although simulations call into question 
what type of gas (inflowing versus IGM) is probed 
along sightlines with metal-poor LLSs. 

This evidence indicates that a statistically significant 
population of low-metallicity LLSs exists at redshift 
z = 3.5-4.5; these absorbers have metallicities consistent 
with being drawn from the IGM and may therefore be an 
observational manifestation of filamentary cold flows pre¬ 
dicted by simulations. Observational and archival pro¬ 
grams that will increase the moderate-resolution sample 
of both the metal-poor and general LLS populations are 
underway. This will allow for more in-depth discussion 
of the metallicity distribution and cosmological impli¬ 
cations, and further coupling to the increasingly more 
detailed analysis of simulated volumes, mapping the dis¬ 
tribution and flow of unprocessed LLS gas relative to 
star-forming galaxies. 
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